%% Uncertainty of Parameter Estimates II

% A Bayesian approach to error estimation

% Estimate the mean from Poisson data and find standard error of the mean


%Create a set of repeat counts from the same experiment
Ntrials=10;
Lambda= 5;
Data=poissrnd(Lambda*ones(Ntrials,1))

%show this data
figure
plot(Data,'o')
xlabel('trials')
ylabel('Counts')

%Calculate the log-likelihood over a range of lambda values
Lvec=linspace(.1,10,100)
LogL=zeros(size(Lvec));

for ll=1:length(Lvec)
   %See  PoissLambdaLogL.m
   [ LL ] = PoissLambdaLogL(Lvec(ll),Data); 
   LogL(ll)=LL; 
end

%Show the LogL plot
figure
plot(Lvec,LogL)
xlabel('Lambda')
ylabel('Neg Log Likelihood')

%Show the likelihood plot
figure
plot(Lvec,exp(-LogL))
xlabel('Lambda')
ylabel('Likelihood')

% Bayes rule:
% P(theta|data)=P(data|theta)*P(theta)/P(data)

% P(data) is integral of P(data|theta)*P(theta) over theta

% By Bayes rule, using a uniform prior, P(theta)=1/L where L is the 
% size of interval of integration. 

L_Lambda=exp(-LogL); %This is P(data|theta)

%Step size for numeric integral
dL=Lvec(2)-Lvec(1);

%This is denominator in Bayes rule.  
NormInt=sum(dL*L_Lambda);

%This is the probability distribution of lambda given our data,
%P(theta|data)
P_LambdaGData=L_Lambda/NormInt;

figure
plot(Lvec,P_LambdaGData)
xlabel('Lambda')
ylabel('P(Lambda|Data)')

% We can approximate this shape using a 1D Gaussian
% the second derivative of a the log of a 1D Gaussian is 1/sigma^2

%The 'fminunc' function works like fminsearch, but also returns the 
% 'Hessian' - the matrix of second derivates.  For a scalar paremeter 
% like we have here, then Hessian is 1/sigma^2

%Setup fminunc and find optimum values. 
Options=optimset();
Theta0=5;
[Theta_Found,~,~,Output,Grad,Hess]=fminunc(@PoissLambdaLogL,Theta0,Options,Data);

%We can give the error on our estimate as the standard error, which is the 
%sigma of our Gaussian distribution
StandardError=sqrt(inv(Hess))

%Plot the distrubution and the Gaussian approximation together
LHess=normpdf(Theta_Found,Lvec,StandardError)
hold on
plot(Lvec,LHess)
legend('P(Lambda|Data)','Gaussian From Hessian')



















